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APPARATUS AND METHOD FOR COMPUTING MULTIPLE INTEGRAL, 
AND RECORDING MEDIUM 
BACKGROUND OF THE INVENTION 
Field of the Invention 

5 The present invention relates to an apparatus and method for computing multiple 
integral, and a recording medium. More particularly, to an apparatus and method for 
high speed multidimensional integration of given integrand functions which are defined 
over a multi-dimensional unbounded domain, and a recording medium in which a 
program for realizing the apparatus and the method for computing multiple integral. 
1 0 Description of the Related Art 

Problems relating to computation of multiple integral on a multi-dimensional 
unbounded domain appear in various technical fields such as risk evaluation in financial 
derivative calculation, evaluating interference noises in multiple-user communication 
systems such as CDMA (Code Division Multiple Access), and the like. 
1 5 In case of numerical integration on 1- or 2-dimensional domain, the domain 
(interval) can be divided into fine sub domains, and an approximated integral can be 
obtained as an average of represented values of an integrand function on each subdomain. 

If the integration domain is in higher dimension, however, the domain must be 
divided into the great number of subdomains. The number of the subdomains must 
20 exponentially increase with respect to the dimensionality of the integration domain to 
maintain the commutative efficiency. This problem is called " curse of dimensionality" . 

Monte Carlo method has been proposed as one of solutions for such the problems. 
Monte Carlo method usually employs pseudo-random numbers (including quasi random 
numbers such as low-discrepancy sequences) on the finite interval which should be 
25 uniformed as possible. 

In a case where integral must be defined in an unbounded domain such as infinite 
domain [-oo, +00], however, it is impossible to normalize uniformed random numbers on a 
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bounded domain after converting the uniformed random numbers on the finite interval 
into uniformed random numbers on the unbounded domain by a linear transformation. 
In this case, a useful non-linear change of variables is required. 

The inverse function method was established by J. von Neumann and S. Ulam in 
5 1947. According to the emthod, it is able to generate random numbers in accordance 
with an arbitrary density function from uniformed random numbers on the finite interval. 
In that year, they carried out computer simulation for tracing fission neutrons, and it was 
successful. Since the method utilizes an inverse function, it is called inverse function 
method. More precisely, the inverse function r\ l (X) is used for generating uniformed 
1 0 random numbers. In this case, an integral function for a density function p is: 
n (X) = Jc p(x) dx (C = {u|-a)<u<X}) 
It has been turned out that random numbers having non-uniformed distribution on a 
bounded domain are available, as chaos theory have been developed recently. Chaotic 
dynamical mapping systems are getting more significant matters for random number 
15 generation. 

Recurrence formula is usually applicable to the random number generation. In this 
case, a rational map is often used as the recurrence formula which generates ideally 
chaotic sequences. The rational map there is a result of the addition theorem of an 
elliptic function (including a trigonometric function). The random number generation 
20 by such the method has the following advantages. 

(1) Non-cyclical random number sequence generation by which the numbers are 
proven to be chaotic (without some exceptions). 

(2) According to sensitivity to initial conditions of chaos, completely different 
random number sequences are available just by slightly changing the seeds (the initial 

25 values to be applied to the recurrence formula). 

(3) Known analytic function acts as the density function expressing random number 
distribution. 
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Known maps which bring the above advantages are: an Ulam-von Neumann map, a 
cubic map, a quintic map, and the like. 

Ulam-von Neumann Map: f(x) = 4x (1-x) 
Cubic Map: f(x) = x (3-4x) 2 
5 Quintic Map: f(x) = x (5-20x+16x 2 ) 2 

Regardless of the different maps above, p(x) = l/«x(l-x)) 1/2 ) represents a density 
function expressing distribution of a random number sequence x[i] (i>0) (x[0], x[l], 
x[2], . . . ) which results from the following recurrence formula: 
x[i+l] = f(x[i]) (i>0) 
1 0 where x[0] =5 (where 0<^<1 ; ^ is an arbitrary initial value) 

Japanese Unexamined Patent Application KOKAI Publication No. HI 0-283344 by 
Umeno et al. discloses a technique for generating random numbers on a bounded domain 
with using those maps, and embodiments where the technique is applied to Monte Carlo 
numerical integration. The following references disclose theoretical backgrounds for the 
1 5 above described application. 

S. M. Ulam and J. von Neumann, Bull. Math. Soc. 53 (1947) pi 120. 
R. L. Adler and T. J. Rivlin, Proc. Am. Math. Soc. 15 (1964) p794. 
K. Umeno, Method of constructing exactly solvable chaos, Phys. Rev. E (1997) vol. 
55, pp. 5280-5284. 

20 The above conventional technique, however, includes the following problems. 

Monte Carlo numerical integration on unbounded domain requires random numbers 
which are distributed in the unbounded domain in accordance with a density function 
given by a known analytic function. Therefore, obtained uniform random numbers must 
be converted into the random numbers distributed in the unbounded domain. To obtain 

25 random numbers distributed in accordance with a desired density function with using the 
aforementioned von Neuman's inverse function method, an integral function of the 
density function must be obtained first, and then, an inverse function of the integral 
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function must be obtained. 

However, it is almost impossible to analytically integrate the density function. 
Moreover, approximation and calculation time accompanying to numerical calculation 
must be considered seriously when obtaining the inverse function of the integral function 
5 resulting from the density function. Since the Monte Carlo method requires the great 
number of random numbers, it is a hard task to calculate the integral with using the Monte 
Carlo method by the inverse function method. 

The numerical integration in a unbounded domain is also required for solving the 
Black-Scholes equation for stock option pricing based on the data of stock prices. If 
1 0 stock price fluctuation follows the Normal distribution (Gaussian distribution), Gaussian 
integration on a unbounded domain is useful for solving the Black-Scholes equation. 

It has been turned out, however, that density function expressing actual distribution 
of stock price fluctuation shows power-law behavior (R. N. Mantegna and H. E. Stanley, 
Nature, vol. 376, pp. 46-49, 1995). Integral of such the density function showing the 
1 5 power-law behavior can not be expressed by elementary functions. Therefore, it is very 
difficult to solve the stock price fluctuation with using the Monte Carlo integration on an 
unbounded domain based on known uniform random numbers. 

Such the problems also appear in evaluation of communication traffic or 
interference noises in CDMA showing power-law behavior in transmission time 
20 distribution or in background noise distribution respectively. 

The application of process to efficient computation of multiple integral on a 
unbounded domain has recently been attracting attention in industry fields. 

It is an object of the present invention to provide an apparatus and a method for 
computing multiple integral, and a recording medium storing a program for realizing the 
25 above, which are suitable for high speed numerical integration of integrands which are 
defined over a multi-dimensional unbounded domain. 

SUMMARY OF THE INVENTION 
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The invention for achieving the above object will now be disclosed in accordance 
with the principal of the present invention. 

An apparatus for computing multiple integral according to a first aspect of the 
present invention computes multiple integral represented by a multidimensional integrand 
5 function A to be integrated with using a vector map f. The vector map f has an 

unbounded support, and converts an m (m>l) -dimensional vector having a real number 
components into an m-dimensional vector having real number components. A 
multidimensional density function p for the limiting distribution resulting from repeatedly 
applying the map f to an m-dimensional vector u is analytically solvable. 
1 0 FIG. 1 is a schematic diagram showing the essential structure of the computing 
apparatus. Components and operations of the computing apparatus will now be 
described with reference to FIG. 1 . 

A computing apparatus comprises a first storage unit, a second storage unit, a first 
computing unit, a second computing unit, an update unit, and an output unit. 
1 5 The first storage unit stores an m-dimensional vector x = (x ls x 2 , . . . , xj. 

The second storage unit stores a scalar value w. 

The first computing unit which computes a vector x' = f(x) = (x',, x' 2 , . . . , x'J 
resulting from applying the vector map f to the vector x being stored in the first storage 
unit. 

20 The second computing unit which computes a scalar value w' = A(x)/p(x) based on 
the vector x being stored in the first storage unit and the scalar value w being stored in the 
second storage unit. 

The update unit updates the value stored in the first storage unit by storing the 
vector x' computed by the first computing unit on the first storage unit, and updates the 
25 value in the second storage unit by adding the scalar value w' computed by the second 
computing unit on the second storage unit. 

The output unit which computes a scalar value s = w/(c+l) based on the scalar value 
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w being stored in the second storage unit when the number of update times by the update 
unit becomes c (c>l), and outputs the scalar value s as a numerical result of the multiple 
integral. 

The present invention is based on the following ergodic property over 
5 multidimensional unbounded domain M: 
J M A(x)dx 

= I M (A(x)/p(x))-p(x)dx 

= (l/N)-E^ N A(XO])/p(X[j]) (N — » co); 

xo+i] = f(xoi) G>0); 

1 0 X[0] is an arbitrary m-dimensional vector 

In this case, computation errors are in the order of (1/N) 1/2 regardless of the number 
of dimensions m even at worst. Here, N is the updating time given by c + 1 . 
According to the latest study by the inventor (K. Umeno "Chaotic Monte Carlo 
Computation: a Dynamical Effect of Random-Number Generations", Japanese Journal of 
1 5 Applied Physics, March 2000), it has been clarified that the computation errors can be in 
the order of 1/N regardless of the number of dimensions m if a condition called Super- 
efficiency is satisfied. Therefore, faster integration is available. 

In the computing apparatus, the scalar value stored in the second storage unit first 
may be a result from dividing a value resulting from applying the function A to the m- 
20 dimensional vector stored in the first storage unit first by a value resulting from applying 
the density function p to the m-dimensional vector stored in the first storage unit first. 

In the computing apparatus, the scalar value stored in the second storage unit first 
may be 0, and 

the output unit may compute a scalar value s' = w/c, and outputs the scalar value s' 
25 as the numerical result of the multiple integral instead of the scalar value s. 

The computing apparatus may further comprise a convergence rate obtainer, a 
vector map selector, and an output controller. 
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The convergence rate obtainer may obtain convergence rate of scalar values 
sequentially output by the output unit while varying the number of update times c for each 
of plural vector maps g l3 g 2 , . . . , g k (k>2) which are prepared as the vector map f. 

The vector map selector may refer to the convergence rates obtained by the 
5 convergence rate obtainer, and select a vector map g h (l<h<k) which shows fastest 
convergence rate. 

The output controller may control the output unit to output the scalar value with 
using the vector map g h as the vector map f and the number of update times c' (c'>c) 
instead of the number of update times c. 
1 0 In the computing apparatus, a multidimensional density function p representing the 
limiting distribution of a vector sequence 

u, f(u), f(f(u)), f(f(f(u))), . . . 
resulting from applying the vector map f to a predetermined m-dimensional vector u = (u l5 
u 2 , . . . , uj for equal to or greater than 0 times, may satisfy the following property of: 
15 P (u) = ri=o m PiOO; 

Pi (v0 ~ c_, kr (I+a) for -co; 
ftOO ~ c_ t |uj" (1+a> for u, -> +oo; 
(a>0, l<i<m, c_>0, c +1 >0) 
In the computing apparatus, the vector map f may be defined as 
20 f(u) = (f^X f 2 (u s ), . . . , UvQ) 

by a function f : (t) = g^d, t)/d, (d>0) which is defined in l<i<m, and the map g, may 
be defined by any one of the following maps % (l<j<8) and a natural number n, (n>2), as 
follows: 

& ((p J (0)) = cp J (n 1 8); 
25 cp.CO) - -sgn(tane) / |tan6| 1/a ; 

cp 2 (6) = -sgn(tanG) x |tanG| 1/a ; 
tp 3 (9) = -sgn(cosO) / Itane| 1/a ; 
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cp 4 (9) = -sgn(cosG) x |tan6| I/a ; 
cp 5 (9) = sgn(cose) / |tan6| 1/a ; 
cp 6 (9) = sgn(cosO) x |tan0| 1/a ; 
q> 7 (9) = sgn(sin6) / |tan0| 1/a ; 
5 cp s (9) = sgn(sine) x |tanG| 1/a ; 
sgn(t) = l fort>0; 
sgn(t) = 0 for t = 0; 
sgn(t) = -1 for t < 0 

The computing apparatus may comprise a convergence rate obtainer, a positive 
1 0 number selector, and an output controller. 

The convergence rate obtainer may define the map f for each of plural positive 
numbers q l5 %. (k>2) prepared as an invariable a, and obtain convergence rates of 

the scalar values sequentially output by the output unit while varying the number of 
update times c. 

1 5 The positive number selector may refer the convergence rates obtained by the 
convergence rate obtainer, and select a positive number q^ (l<h<k) which shows the 
fastest convergence rate. 

The output controller may define the map f with using the positive number q^ as the 
invariable a, and control the output unit to output the scalar values with using the number 
20 of update times c' (c'>c) instead of the number of update times c. 

The computing apparatus may comprise a convergence rate obtainer, a map selector, 
and an output controller. 

The convergence rate obtainer may define the map g t with using plural ones of the 
maps cp j3 and obtain convergence rates of the scalar values sequentially output by the 
25 output unit while varying the number of update times c. 

The map selector may refer to the convergence rates obtained by the convergence 
rate obtainer, and select one of the maps cpj which shows the fastest convergence rate. 
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The output controller may define the map & with using the map cpj selected by the 
map selector, and control the output unit to output the scalar values with using the number 
of update times c' (c'>c) instead of the number of update times c. 

The computing apparatus may comprise a convergence rate obtainer, a natural 
5 number selector, and an output controller. 

The convergence rate obtainer may define the map g, relating to each of plural 
natural numbers p b p 2 , . . . , p k (k>2) as the natural numbers % and obtain convergence 
rates of the scalar values sequentially output by the output unit while varying the number 
of update times c. 

1 0 The natural number selector may refer to the convergence rates obtained by the 
convergence rate obtainer, and select a natural number p h (l<h<k) which shows the fastest 
convergence rate. 

The output controller may define the natural number map g, with using the natural 
number p h as the natural number n i; and control the output unit to output the scalar values 
1 5 with using the number of update times c' (c'>c) instead of the number of update times c. 
The output unit of the computing apparatus may compute the scalar value s each 
time the update unit carries out update, compare the latest scalar value s with the former 
scalar value which is computed at former update, and output the latest scalar value s if a 
result of the comparison satisfies a predetermined condition for terminating the 
20 computation. 

A method for computing multiple integral according to a second aspect of the 
present invention computes multiple integral of a multidimensional integrand function A 
to be integrated with using a vector map f, a first storage unit which stores an m- 
dimensional vector x = (x l9 x 2 , . . . , xj, and a second storage unit which stores a scalar 
25 value w. 

The vector map f has an unbounded support, and converts m (m>l) -dimensional 
vector having real number components into m-dimensional vector having real number 



10 

components. A multidimensional density function p for the limiting distribution 
resulting from repeatedly applying the map f to m-dimensional vector u is analytically 
solvable. 

The computing method comprises the first computing step, the second computing 
5 step, the updating step, and the outputting step. 

The first computing step computes a vector x" = f(x) = (x' la x' 2 , . . . , x' J resulting 
from applying the vector map f to the vector x being stored in the first storage unit. 

The second computing step computes a scalar value w' = A(x)/p(x) based on the 
vector x being stored in the first storage unit and the scalar value w being stored in the 
1 0 second storage unit. 

The updating step updates the value stored in the first storage unit by storing the 
vector x' computed by the first computing unit on the first storage unit, and updating the 
value in the second storage unit by adding the scalar value w' computed by the second 
computing unit to a value to be stored on the second storage unit. 
1 5 The outputting step computes a scalar value s = w/(c+l) based on the scalar value w 
being stored in the second storage unit when the number of update times by the update 
unit becomes c (c>l), and outputs the scalar value s as a numerial result of the multiple 
integral. 

In the computing method, the scalar value stored in the second storage unit first may 
20 be a result from dividing a value resulting from applying the function A to the m- 

dimensional vector stored in the first storage unit first by a value resulting from applying 
the density function p to the m-dimensional vector stored in the first storage unit first. 

In the computing method, the scalar value stored in the second storage unit first may 
be 0, and 

25 the outputting step may compute a scalar value s' = w/c, and output the scalar value 
s' as the numerical result of the multiple integral instead of the scalar value s. 

The computing method may further comprises the convergence rate obtaining step, 
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and the vector map selecting step. 

The convergence obtaining step may obtain convergence rate of scalar values 
sequentially output by the outputting step while varying the number of update times c for 
each of plural vector maps g ls g 2 , . . . , g k (k>2) which are prepared as the vector map f. 
5 The vector map selecting step may refer to the convergence rates obtained by the 
convergence rate obtainer, and select a vector map g h (l<h<k) which shows fastest 
convergence rate. 

A vector map g h may be used as the vector map f. 

The outputting step may output the scalar value with using the vector map g h as the 
1 0 vector map f and the number of update times c 1 (c'>c) instead of the number of update 
times c. 

In the computing method, a multidimensional density function p representing the 
limiting distribution of a vector sequence 

u, f(u), f(f(u)), f(f(f(u))), . . . 
1 5 resulting from applying the vector map f to a predetermined m-dimensional vector u = (u ls 
u 2 , . . . , uj for equal to or greater than 0 times, may satisfy the following property of: 

p(u) = n 1 =o m P i (u 1 ); 

Pi(Ui) ~ c_, luj^ for u, -co; 
pXuJ ~ c_, | Ui r (1+a) for u, -* +co; 
20 (a>0, l<i<m, c_>0, c + >0) 

In the computing method, the vector map f may be defined as 

f(u) = (f^uO, f 2 (u s ), . . . , fjuj) 
by a function f,(t) = g,(d, t)/dj (d >0) which is defined in l<i<m, and the map g, is 
defined by any one of the following maps cp 3 (l<j<8) and a natural number n t (n>2), as 
25 follows: 

g l (cp J (e)) = (p j (n i e); 

(p^G) = -sgn(tan6) / |tan6| 1/a ; 
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(p 2 (9) = -sgn(tan0) x |tanO| 1/a ; 

<p 3 (9) = -sgn(cosG) / |tan9| 1/a ; 

<p 4 (9) = -sgn(cosG) x |tan9| 1/a ; 

<p 5 (9) = sgn(cos6) / |tan0| 1/a ; 
5 (p 6 (0) = sgn(cos9) x |tan9| 1/a ; 

(p 7 (0) = sgn(sinG) / |tan9| 1/a ; 

<p 8 (0) = sgn(sin9) x |tan9| 1/a ; 

sgn(t) = 1 for t > 0; 

sgn(t) - 0 for t = 0; 
10 sgn(t) = -1 for t < 0 

The computing method may further comprise the convergence rate obtaining step, 
the positive number selecting step, and the output controlling step. 

The convergence rate obtaining step may define the map f for each of plural positive 
numbers q l9 . . . , % (k>2) prepared as an invariable a, and obtain convergence rates of 
1 5 the scalar values sequentially output by the output unit while varying the number of 
update times c. 

The positive number selecting step may refer to the obtained convergence rates, and 
select a positive number q h (l<h<k) which shows the fastest convergence rate. 

The output controlling step may define the map f with using the positive number 
20 as the invariable a, and control output of the scalar values with using the number of 
update times c' (c'>c) instead of the number of update times c. 

The computing method may further comprises the convergence rate obtaining step, 
the map selecting step, and the output control step. 

The convergence rate obtaining step may define the map g, with using plural ones of 
25 the maps q> and obtain convergence rates of the scalar values sequentially output by the 
output step while varying the number of update times c. 

The map selecting step may refer to the obtained convergence rates, and select one 
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of the maps qx, which shows the fastest convergence rate. 

The output controlling step may define the map g { with using the selected map cp ; 
selected, and control output of the scalar values with using the number of update times c' 
(c'>c) instead of the number of update times c. 
5 The computing method may further comprises the convergence rate obtaining step, 
the natural number selecting step, and the output controlling step. 

The convergence rate obtaining step may define the map g, relating to each of plural 
natural numbers p l9 p 2 , . . . , p k (k>2) as the natural numbers n i3 and obtain convergence 
rates of the scalar values sequentially output by the outputting step while varying the 
1 0 number of update times c. 

The natural number selecting step may refer to the obtained convergence rates, and 
select a natural number p h (l<h<k) which shows the fastest convergence rate. 

The output control step may define the natural number map ^with using the natural 
number p h as the natural number n,, and control output of the scalar values with using the 
1 5 number of update times c' (c'>c) instead of the number of update times c. 

In the computing method, the outputting step may compute the scalar value s each 
time the updating step carries out update, compares the latest scalar value s with the 
former scalar value which is computed at former update, and output the latest scalar value 
s if a result of the comparison satisfies a predetermined condition for terminating the 
20 computation. 

A program for realizing the apparatus and method for computing multiple integral 
according to the present invention may be stored in a recording medium such as a 
compact disc, a floppy disk, a hard disk, a magneto-optical disk, a digital versatile disc, a 
magnetic tape, or a semiconductor memory. 
25 The program stored in the recording medium according to the present invention may 
be executed by a data processor being equipped with a storage device, a computing unit, 
output device, and the like, such as a general purpose computer, a video game device, a 
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PDA (Personal Data Assistants), and a cellular phone, to realize the above described 
apparatus and method for computing multiple integral. 

The recording medium storing the program according to the present invention may 
be distributed or merchandized being separated from the data processor., 
5 BRIEF DESCRIPTION OF THE DRAWINGS 

These objects and other objects and advantages of the present invention will become 
more apparent upon reading of the following detailed description and the accompanying 
drawings in which: 

FIG. 1 is a schematic diagram showing the essential structure of an apparatus for 
1 0 computing multiple integral according to the present invention; 

FIG. 2 is a schematic diagram showing the essential structure of a data processor in 
which the apparatus for computing multiple integral according to the present invention is 
installed; 

FIG. 3 is a flowchart showing the steps of the method of computing multiple 
1 5 integral according to the present invention; 

FIG. 4 is a graph exemplifying a map £ which defines a vector map f; 
FIG. 5 is a graph showing a density function of the limiting distribution in a scalar 
sequence obtained by applying the map f, repeatedly; 

FIG. 6 is a graph showing means square of errors in computing steps where 1- 
20 dimensional integrand function to be numerically integrated in infinite domain (-oo, +co); 

FIG. 7 is a graph showing mean square deviation of errors in computing steps where 
2-dimensional function to be integrate is numerically integrated in unbounded domain 

(—oo 3 +oo) x (— oo 3 +oo); 

FIG. 8 is a flowchart showing steps of the method of computing multiple integral 
25 according to the present invention; 

FIG. 9 is a flowchart showing steps of selecting a vector map according to the 
present invention; and 
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FIG. 10 is a flowchart showing steps of selecting a vector map according to the 
present invention. 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
Embodiments of the present invention will now be described. One skilled in the 
5 art may be able to propose modified embodiments each of which include all or some 
elements described in the following embodiments of the present invention. Such the 
modified embodiments will be included in the scope of the present invention, because the 
following embodiments of the present invention do not limit the scope of the present 
invention but just explain the present invention. 

1 0 First Embodiment 

FIG. 2 is a schematic diagram showing the essential structure of a data processor 
such as a general purpose computer which realizes the apparatus for computing multiple 
integral according to the present invention. A first embodiment of the present invention 
will now be described with reference to FIG. 2. 

1 5 A data processor 30 1 is controlled by a CPU (Central Processing Unit) 302. A 
main storage 303 such as RAM (Random Access Memory) temporarily holds data. An 
external storage 304 such as a hard disk, a floppy disk, a CD-ROM (Compact Disc Read 
Only Memory), a magnetic tape medium, and a magneto-optical disk stores programs to 
be executed by the CPU 302. 

20 When the data processor 30 1 is booted, the CPU 302 executes an initial program 
loader (IPL) from ROM (Read Only Memory) 308. According to the execution of the 
IPL, programs relating to the operating system, application programs, and/or the like are 
loaded to the main storage 303 from the external storage 304 to be executed. 

Results from execution of such the programs are stored on the external storage 304 

25 as files or displayed on a display unit 305 such as a CRT (Cathode Ray Tube) display and 
a LCD (Liquid Crystal Display). A user of the data processor 301 operates input devices 
306 such as a mouse and a keyboard to input instructions to the data processor 301. 
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The main storage 303 acts as the first storage unit and the second storage unit. 
The CPU 302 is equipped with ALU (Arithmetric Logic Unit) and a coprocessor for 
numerical computation (not show) which act as the first and second computing unit. 
The CPU 302 collaborates with the main storage 303 to act as the update unit and 
5 the output unit. 

The external storage 304 acts as the recording medium storing a program according 
to the present invention. 

Since those functions are realized by numerical computation and memory update as 
described below, the structure of the present invention can be realized by an electronic 
1 0 circuit designed for the above functions instead of the general purpose computer. The 
scope of the present invention includes such the modification. 
Computing Process 

FIG. 3 is a flowchart showing steps to be executed by the computing apparatus 
according to the present invention, that is, showing process flow of the steps in the 
1 5 computing method according to a first embodiment of the present invention. 

The main storage 303 (typically, RAM) previously has an area for m-dimensional 
variable x (first storage unit), an area for scalar variable w (second storage unit), another 
areas for m-dimensional variable x' and for scalar variable w' which are prepared for 
temporary storage areas, an area for number of update times c, and an area for integration 
20 result variable s. Those areas will be denoted simply as variable x, variable w, variable 
x', variable w', variable c, and variable s respectively. And values to be stored in those 
areas will be denoted as x, w, x', w', c, and s respectively. 

"A" denotes multidimensional integrand function to be integrated over an 
unbounded domain M = (-co, +oo) x . . . x (-co, +oo), "f ' denotes vector chaotic map 
25 satisfying later-described conditions, and "p" denotes a density function of the the 
limiting distribution after applying the vector chaotic map f repeatedly. 

The steps will now be described in detail with reference to FIG. 3. First, the CPU 
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302 stores a value 0 on the variable c (step S301). 

Then the CPU 302 stores an arbitrary m-dimensional initial vector on the variable x 
(step S302). 

The CPU 302 computes a scalar value A(x)/p(x), and stores the result on the 
5 variable w (step S303). 

The CPU 302 computes a vector value f(x), and stores the result on the variable x' 
(step S304). 

The CPU 302 calculates the scalar value A(x)/p(x), and stores the result on the 
variable w' (step S3 05). 
1 0 Then, the CPU 302 updates the variable x by storing x' on the variable x (step S306). 
The CPU 302 then updates the variable w by adding w' to the variable w (step 

5307) . 

The CPU 302 further updates the variable c by adding 1 to the variable c (step 

5308) . 

1 5 The CPU 302 determines whether c exceeds predetermined number of repeat times 
C or not (step S309). If c does not exceed C yet (step S309: No), the flow returns to step 
S3 04. 

On the contrary, if c exceeds the predetermined number of repeat times C (step 
S309: Yes), the CPU 302 computes a scalar value w/(c+l), and stores the result on the 
20 variable s (step S3 10). And the process flow is terminated, "s" represents the 
numerical result of integration. 
Examples of Vector Maps 

A vector map f being applicable to the present invention is: 
f(u) = (f 1 (u 1 ),f 2 (u 2 ),...,f m (u m )) 
25 where, m-dimensional vector u = (u x , u 2 , . . . , uj, and 

• map f t> f 2 , . . . , 4 
FIG. 4 exemplifies maps f- (l<i<m) which define a vector map f. 
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Each map £ is called as chaotic mapping. As aforementioned, those are available 
from n-multiplication formulas of a tangent function (n>2). The inventor wrote the 
mathematical theorems and demonstrations for the above in: 

Ken Umeno, Superposition of chaotic processes with convergence to Levy's stable 
5 law, Physical Review E, vol. 58 No. 2. The American Physical Society, Aug. 1998. 

F2(x) in FIG. 4 represents the map £ when parameters are n=2, d=l, and a=l, and q> 
is used. 

F2*(x) in FIG. 4 represents the map £ when parameters are n=2, d=l, and a=l, and 
cp 5 is used. 

1 0 F3(x) in FIG. 4 represents the map fj when parameters are n=3, d=l, and a=l, and cp 
is used. 

F3*(x) in FIG. 4 represents the map £ when parameters are n=3, d=l, and a=l, and 
cp 5 is used. 

Explicit forms of those functions are available as follows. 
15 F2(x) = 2|xj/(l-x 2 ) 

F2*(x) = l/2(lx|-l/|x|) 
F3(x) = x|x 2 -3|/(l-3x 2 ) 
F3*(x) = -x(x 2 -3)/|(l-3x 2 )| 

FIG. 5 shows densities obtained by a density function representing the limiting 
20 distribution of a scalar sequence t, £(t), £(£(t)), £($(£(1:))), . . . which is obtained by 
applying the map £ repeatedly. 

FIG. 5 shows plotted density where a is varied while d is fixed to 1 in the density 
function (a represents a in FIG. 5). As mentioned in the above described writing, the 
density function is analytically solvable as: 
25 p(x) = ad a |x| a -'/(7i(l+d 2a |x| 2a )) 

It is obvious from FIG. 5 that shape of the density function p radically changes as a 
varies, regardless of n and q> } . 
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Experiment Result 1 

FIG. 6 is a graph showing mean square deviation of computation error over the 
number of computing steps N = c + 1 (c is the number of update times), hi this case, the 
computation is to integrate 2 different 1 -dimensional integrand functions (scalar maps) 
5 A(-), A'(-) to be integrated with using various chaotic maps f, (P3(X)-P1(X) and P3*(X)- 
P 1 *(X) represent A(x) and A'(x) respectively in FIG. 6). The 1 -dimensional integrand 
functions are expressed as follows: 

A(x) = -4x 2 sgn(x)/(7<l+x 2 ) 5/2 ); 
A(x) = -4|x|sgn(x)/«l+x 2 f 2 ) 
1 0 An experiment where the above described F3(x) and F3*(x) are applied to those 
integrand functions was carried out to obtain average errors after independent trials were 
carried out 1,000 times for each step number N. 

Since it has been analytically known that the result will be 0 after integrating A(-) 
and A'(0 on the infinite domain (-co, +co), obtained result in the numerical integration 
1 5 will be the computation error. 

Mean square deviation of errors radically decreases in 1/N 2 order where N is the 
number of the steps, when F3(x) is applied to the integrand function A(-)p(0 and F3*(x) is 
applied to another integrand function A'(-)pC) respectively. 

On the other hand, decrease rate of the square deviation of errors is loosened when 
20 F3*(x) is applied to the integrand function A(-)p(0 and F3 (x) is applied to another 

integrand function A'Qp(') respectively, that is, the means square deviation decreases in 
1/N order where N is the number of the steps. 

In conclusion, faster convergence than ordinary Monte Carlo computation is 
available when the most suitable chaotic map is selected as a random number generator. 
25 Experiment Result 2 

FIG. 7 is a graph showing mean square deviation of computation error over the 
number of computing steps N. In this case, the computation is to integrate 2- 
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dimensional integrand function A(-, ■) on 2-dimensional unbounded domain M = (-00, 
+00) x(-oo, +00) by the technique of the present invention. The 2-dimensional function is 
expressed as follows: 
A(x, y) 

5 = ((l-x 2 )/(l+x 2 )) ((l-Sy^/Cl+y 2 ) 3 ' 2 ) sgn(y) - sgn(x) (l+x 2 ) 1/2 |y| sgn(y) (l+y 2 ) 1 ' 2 

FIG. 7 shows results of four cases where ft, f 2 ) = (F2, F2), (F2, F3), (F3, F2) or (F3, 
F3) . As well as the experiment 1, average computation errors after carrying out 
independent trials 1,000 times for each step number N are obtained. 

Since it has been analytically known that J M A(x, y)dxdy = 0, the result of numerical 
1 0 integration will be the computation error as well as the above indicated case. 

According to FIG. 7, of four cases, the case where F2(x) is selected as the map fi 
and F3(x) is selected as the map f 2 shows rapid decrease of the mean square deviation 
over the step numbers N. In this case the decrease rate is in 1/N 2 order. Other cases 
show gentler decrease of the mean square deviation in the errors. In this case, the 
1 5 decrease rate is in 1/N order. 

It is also obvious from the experiment result 2 that faster convergence of 
computation errors than ordinary Monte Carlo method is available by selecting a suitable 
chaotic map as the random number generator from various chaotic maps. 

Instead of the chaotic map f described above, the following definition may be 
20 applicable to this embodiment. 

^u) = (f 1 (u 1 ,u 2 ,...,uj s 

f 2 (u ls u 2 , . . . , uj, . . . , 
4( Ul , u 2 , . . . , uj); 
u = (u 1; u 2 , . . . , uj 

25 The present invention is able to employ the above definition if random variable 
density in the limiting distribution of vector sequence resulting from applying f to 
predetermined vector repeatedly, has an explicit form. 
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In addition to the above described Ulam-von Neumann map, cubic map, and quintic 
map, Chebyshev map or generalized Chebyshev map is also applicable as f r 
Second Embodiment 
FIG. 8 is a flowchart showing steps to be executed by the computing apparatus 
5 according to the present invention, that is, showing process flow of the steps in the 
computing method according to a second embodiment of the present invention. 
Variables, maps, functions, and the like in the second embodiment are the same as those 
in the first embodiment. 

The process flow will now be described with reference to FIG. 8. First, the CPU 
1 0 302 stores a value 0 on the variable c (step S801). 

Then the CPU 302 stores arbitrary m-dimensional initial vector on the variable x 
(step S802). 

The CPU 302 stores 0 on the variable w (step S803). 

The CPU 302 computes a scalar value A(x)/p(x), and updates the variable w by 
1 5 adding the resultant scalar value to the variable w (step S804). 

The CPU 302 further updates the variable c by adding 1 to the variable c (step 
S805). 

Then, the CPU 302 determines whether the variable c exceeds the number of repeat 

times C or not (step S806). 
20 If the variable c does not exceed the number of repeat times C yet (step S806: No), 

the CPU 302 computes m-dimensional vector value f(x) with using the variable x and a 

vector map f, and stores the result on the variable x to update it (step S807), then the flow 

returns to step S804. 

On the contrary, if the variable c exceeds the predetermined the number of repeat 
25 times C (step S806: Yes), the CPU 302 computes a scalar value w/c, and stores the result 

on the variable s (step S809). Then, the process flow is terminated. In this case, the 

variable s represents the result of integration. 
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The second embodiment features that required areas for storing variables are less 
than those required in the first embodiment. 

Third Embodiment 

As described above, if the chaotic map f shows a property called "Super-efficiency", 
5 faster error convergence is available by Monte Carlo computation. However, some 
integral function A have difficulties in selecting a suitable vector map f to show Super- 
efficiency. A third embodiment features automatically selection of the vector map. 
FIG. 9 is a flowchart showing process flow for selecting steps of the vector map f. 
Plural vector maps to be selected for numerical integration are prepared (step S901). 
1 0 For example, various maps such as the aforementioned F 2 , F 2 *, F 3 , and F 3 * having 

different parameters are prepared as maps for computing components of a vector map f. 

The above described method of multiple integral computing is applied to each of the 
prepared vector maps for predetermined times cl and c2 (cl<c2), to obtain approximates 
of integration(step S902). This step is carried out for previously obtaining convergence 
1 5 by each map. 

The obtained plural values resulting from the numerical integration are compared 
with each other to select a vector map which shows fastest convergence (step S903). 
The convergence rate may be determined by obtaining differences between the values 
after the numerical integration by computing c a times and another values after the 

20 numerical integration by computing c 2 times, vector map by vector map. Of the vector 
maps, one having the least absolute value of the obtained difference, which means that it 
shows faster convergence, is selected. 

Or, computing c 2 times may be omitted. In this case, a vector map which outputs a 
value closes to average of the obtained values after the numerical integration, should be 

25 selected. 

Finally, the above described multiple integral computation is carried out over the 
selected vector map for C (c 2 <C) times to obtain the numerical integration (step S904). 
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According to the third embodiment, a vector map which shows faster convergence is 
automatically selected even if it is unknown which vector map shows Super-efficiency for 
faster convergence, thus fast numerical integration is realized. 

Fourth Embodiment 

5 Monte Carlo method often requires predetermined number of random numbers 
(number of repeat times, number of update times according to the present invention) 
before computing numerical integration. The following fourth embodiment features that 
repeat of computation is stopped when the error becomes smaller than a predetermined 
value. 

1 0 For example, a scalar value w/c is computed at each computation (first time, second 
time, . . . ) executed repeatedly as described in the second embodiment. Hereinafter, the 
resultant values will be represented as: 

Sj, s 2 , . . . , s,., . . . 

In the above described embodiments, the computation is terminated when the 
1 5 number of repeat times exceeds the predetermined number of update times C. In this 
embodiment, predetermined error rate s and least number of repeat times L (L>2) are 
determined. And, a scalar value s, will be output as a result of numerical integration if 
the following condition for termination is satisfied: 
i > L; 

20 ls.-s._J < s |sj 

FIG. 10 is a flowchart showing steps of computing numerical integration according 
to the fourth embodiment. Variables, maps, and functions in this embodiment are the 
same as those used in the above embodiments. Additionally, an area for storing a scalar 
variable s' is prepared in the main storage 303. 
25 The process flow will now be described with reference to FIG. 10. The CPU 302 
stores a value 0 on the variable c (step S1001). 

Then the CPU 302 stores arbitrary m-dimensional initial vector on the variable x 
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(step S1002). 

The CPU 302 further stores 0 on the variable w (step S1003). 
Then, the CPU 302 computes a scalar value A(x)/p(x) with using the value stored in 
the variable x, and adds the result to the variable w, thus the variable w is updated (step 
5 S1004). 

The CPU 302 adds 1 to the variable c to update it (step S1005). 
Then the CPU 302 computes a scalar value w/c, and stores the result on the variable 
s' (step SI 006). 

The CPU 302 determines whether the termination condition "c>L n |s'-s| < |es'|" is 
1 0 satisfied or not (step S 1007). 

If the condition is not satisfied (step S1007: No), the CPU 302 stores s' on the 
variable s to update it (step SI 008), then, computes m-dimensional vector value f(x) with 
using the variable x and the vector map f, and stores the result on the variable x to update 
it (step S 1 009). The flow returns to step S 1 004. 
15 On the contrary, if the condition is satisfied (step SI 007: Yes), the CPU 302 stores 
the value s' on the variable s (step S 1 0 1 0), and the process flow is terminated. In this 
case, the value s represents the result of the integration. 

Arbitrary desired value may be applicable to the error rate s, for example, 0.01 (1%) 
or 0.001 (0.1%). 

20 The least number of repeat times L may depend on available computing time or type 
of the integral function A. The least number of repeat times L may be, for example, in a 
range of 100 to 1,000 times. 

Or, the following example may be applicable to this embodiment. If the 

f 

termination condition is satisfied successively k times (for example, 2 to 50 times) in the 
25 computing loop, that is: 

i + k > L; 

js 1+ - Si | < s [s i+1 j; 
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l S i+2 S i+ll — S l S i+2^ 
|S 1+k -S 1+k _!| < 8 |S i+k |, 

the repeat is terminated, and s i+k is output as the result of the numerical integration. 
5 Such the modification may be included in the scope of the present invention. 

Fifth Embodiment 

In the above embodiments, unbounded domain is represented as (-co, +co) x . . . x 
(-oo, +co), and integral is obtained over whole domain. The following fifth embodiment 
discloses a technique for computing numerical integration over unbounded domain M 
1 0 other than the above. Such the unbounded domain M may be, for example, (0, +oo); 
(0, +oo) x (0, +oo); or (0, 1) x (-co, +ao). 

To realize the computation over such the unbounded domain M, it may determine 
whether the updated vale of the variable x is included in the unbounded domain M or not, 
after step S807 described in the second embodiment. If the updated value is included in 
1 5 there, the flow may return to step S804. On the contrary, if the updated value is not 
included there, the flow may returns to step S805. 

According to this modification, fast numerical integration over arbitrary unbounded 
domain M is also available as well as the above indicated embodiments. 

According to the present invention as described above, it is able to provide an 
20 apparatus and a method for computing multiple integral, and a recording medium storing 
a program for realizing the apparatus and the method. 

Those are applicable to various industrial fields such as risk evaluation in financial 
engineering, interference noise evaluation in CDMA technology for mobile 
communications or optical telecommunications technique, and traffic evaluation in 
25 telecommunications such as Internet. 

Various embodiments and changes may be made thereunto without departing from 
the broad spirit and scope of the invention. The above-described embodiments are 
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intended to illustrate the present invention, not to limit the scope of the present invention. 
The scope of the present invention is shown by the attached claims rather than the 
embodiments. Various modifications made within the meaning of an equivalent of the 
claims of the invention and within the claims are to be regarded to be in the scope of the 
5 present invention. 

This application is based on Japanese Patent Application No. HI 1-370277 filed on 
December 27, 1999 and including specification, claims, drawings and summary. The 
disclosure of the above Japanese Patent Application is incorporated herein by reference in 
its entirety. 



